#include <fstream>
#include <istream>
#include <sstream>
#include <vector>
#include <deque>
#include <string>
#include <iostream>



#define DEBUG_LEVEL 6
#include <debug_func.h>

#include <seq_pairs.h>
#include <blast_results.h>
#include <bcpp_read_seq.h>
#include <codeml_results.h>
#include <progress_indicator.h>
#define _MY_STL_MACROS
#include <stlhelper.h>
#include <ultoa.h>
#include <print_error.h>

#include <boost/foreach.hpp>

#include "get_arg.h"


using std::ifstream;
using std::istream;
using std::istringstream;
using std::string;
using std::cerr;
using std::cin;
using std::cout;
using std::runtime_error;
using BioCPP::t_collect_sequences;
using BioCPP::t_collect_sequences_report_errors_to_strm;



//________________________________________________________________________________________

//		do_parse_kaks

//________________________________________________________________________________________
struct do_parse_kaks
{
	t_collect_sequences&	seq_db;
	float					dS;
	float					coverage;
	float					percent_id;
	do_parse_kaks(t_collect_sequences&	sequences_, float ds_, 
				  float coverage_, float percent_id_):
				seq_db(sequences_), dS(ds_), coverage(coverage_), percent_id(percent_id_)
				{	}
	void operator()(istream& is, unsigned)
	{
		// get data for one protein pair
		string acc1, acc2;
		is >> acc1 >> acc2;
		
		t_codeml_results dndsdata;
		is >> dndsdata;
		t_blast_results aligndata;
		is >> aligndata;

		if (!is)
		{
			if (!is.eof())
				throw runtime_error("The kaks distances data is not in the correct format");
			return;
		}

		// make sure > coverage
		if (float(aligndata.aligned_residues) / (std::max(aligndata.len1, aligndata.len2))
			< coverage)
			return;

		// make sure > dS
		if (dndsdata.dS < dS)
			return;
		
		// get protein and gene indices
		unsigned seq_id1, seq_id2;
		if (!seq_db.get_index_of_name(acc1, seq_id1))
			throw runtime_error("The sequence for ["  + acc1 + "] is missing.");
		if (!seq_db.get_index_of_name(acc2, seq_id2))
			throw runtime_error("The sequence for ["  + acc2 + "] is missing.");

		string seq1 = seq_db.sequence(seq_id1);
		string seq2 = seq_db.sequence(seq_id2);
		insert_gaps_into_sequence_pair(aligndata, seq1, seq2,
										false, acc1, acc2);
		unsigned identical = 0;
		for (unsigned i = 0; i < seq1.size(); ++i)
			if (seq1[i] == seq2[i])
				++identical;
		aligndata.percent_id = identical * 1.0 / seq1.length();

		if (aligndata.percent_id < percent_id)
			return;
		
		cout
			<< acc1 << "\t" << acc2 << "\t" 
			<< dndsdata << "\t"
			<< aligndata << "\n";
	}
};






//________________________________________________________________________________________

//		read_sequences

//________________________________________________________________________________________
void read_sequences(t_collect_sequences& seq_db, t_program_args& args)
{
	{
		ifstream f_sequences(args.file_seq.c_str());
		if (!f_sequences)
			throw runtime_error("Could not open " + args.file_seq +
										" to retrieve cDNA sequences");

		args.VERBOSE()	<< "\tReading protein sequences  from [ "
						<< args.file_seq << " ]...\n";
		t_progress_indicator dots(args.VERBOSE(), 1000);
		t_collect_sequences_report_errors_to_strm report_errors(args.CERR());
		seq_db.init_from_stream(f_sequences, ".*", dots, report_errors);
	}
	args.VERBOSE() << "\tDone...\n";
	args.VERBOSE() << "\t" << seq_db.size() << " sequences\n";
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_program_args args;

	try
	{

		debug_func_cerr(5);

		// don't use C streams
		std::ios::sync_with_stdio(false);


		// get command line arguments
		if (!process_command_line_options(argc, argv, args))
			return 1;


		//
		// read sequences
		//
		t_collect_sequences seq_db;
		read_sequences(seq_db, args);

		do_parse_kaks parse_kaks(seq_db, args.dS,  args.coverage,  args.percent_id);
		seq_pairs::read_data(parse_kaks, args.file_input, args.VERBOSE());
	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}

